
# Climate change impact projections for DJO 2012

# Trying our best to reproduce their impact recipe in Appendix III of their 2008 NBER working paper.  Evaluating over bootstrapped coefficients from specification in Table 2, Column 3. They appear to ignore any precip effects so we will too. 

rm(list=ls())  #clear working memory
wd <- "C:\\Users\\Elisa Cascardi\\Desktop\\BDLMS_Replication\\"  #insert into quotation marks the working directory to location where replication file was unzipped. E.g /Documents/ReStat_Replication/ 	Make sure to include the final backslash. You may need to replace single backslash with double backslash depending on operating system.

"%&%"<-function(x,y)paste(x,y,sep="")  #define a function for easy string pasting
setwd(wd%&%"data\\DJO\\")  #navigate to DJO data directory


# Bring in the data: boostrapped regression coefficients created in "DJO_uncertainty.do", and country-level population-weighted climate change projections
boot=read.csv("boot_djo.csv")  #boostrapped reg coefficients from DJO specification in Table 2, Col 3. These are generated in the stata script
param=read.csv("djo_paramdata.csv") #parameters needed to calculate climate growth impacts, following DJO method in 2008 NBER paper. These are generated in the stata script
temp=read.csv("DJO_tempchg_popweight.csv")  #temperature projections

#collect scenario and model names
nms=strsplit(names(temp)[3:104],"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][1]) #scenario
	mod=c(mod,nms[[i]][2]) #climate model
	per=c(per,nms[[i]][3]) #mid- or end-of-century
}


# DJO's future country-level growth rates depend on a convergence parameter, which relates to the difference between a country's log gdp in a given year and the US's gdp in a given year. This is implemented as commented below. 
# Equation for country-specific growth rate in a given year is given in Appendix 3 of DJO 2008 NBER paper
#	g_it = base_us_growth + beta_i (ln US_t - ln_cty_t) + (temp_beta)*tchg*(country==poor)
# We interate year by year for each climate model (tchg) and each temp_effect estimate from the bootstrap.

# Set initial parameters
param$beta[param$fips=="US"]=0
param=param[is.na(param$rgdpl)==F&is.na(param$beta)==F,]
loc=match(param$fips,temp$fips)  #location of DJO countries in climate change projections
param=param[is.na(loc)==F,]  #we don't match to UAE for some reason, but it is not in poor country sample
loc=loc[is.na(loc)==F]
beta=param$beta
beta[beta<0]=0  #following their instructions, set this to zero if it's negative
base=rep(mean(param$g_us),length(beta))  #base rate is US growth if growing faster than US
base[param$beta<0]=param$g[param$beta<0] #otherwise base rate is that country's rate of growth
inipoor=param$rgdpl<3170  #countries that were initially poor

# Initialize parameters and run a no climate change scenario
lus_gdp=log(param$rgdpl[param$fips=="US"])  #initialize US gdp
gdp=param$rgdpl
lgdp=log(gdp)
poor=(gdp<3170)  #redefine every period

#first run the no climate change baseline. i.e. tchg=0
for (i in 4:90) {  #DJO sample is through 2003, so starting at 2004 and running through 2090 (i.e. 2080-2100)
	git = base + beta*(lus_gdp - lgdp)  #growth rate for that year
	gdp= gdp*(1+git/100)
	lus_gdp=log(gdp[param$fips=="US"])  #re-initialize US gdp
	lgdp=log(gdp)
	poor=(gdp<3170)  #redefine who's poor
	if (i==50) {mgdp_nocc_mid=mean(gdp[inipoor])}  #saving a mid-century estimate of median per cap gdp in initially poor countries
}
mgdp_nocc_end=mean(gdp[inipoor])  #save end-century estimate


# Now loop over bootstrapped coefficients, and over climate change scenarios
impact=array(dim=c(dim(boot)[1],102))  #102 is total number of models*scenarios 

for (b in 1:dim(impact)[1]) {  #looping over bootstraps
	tbeta=boot[b,2] #use the point estimate for now	
for (tm in 1:dim(impact)[2]) {
	tchg=temp[loc,tm+2]  # first two columns are names
	#reinitialize parameters
	lus_gdp=log(param$rgdpl[param$fips=="US"])  #initialize US gdp
	gdp=param$rgdpl
	lgdp=log(gdp)
	poor=(gdp<3170) 
	if (per[tm]=="mid") {yrs=4:50} else {yrs=4:90}
	for (i in yrs) {   
		git = base + beta*(lus_gdp - lgdp) + tbeta*(tchg*i/max(yrs))*poor  #growth rate for that year. tchg is calculated such that it is relative to base year growth rate.
		gdp= gdp*(1+git/100)
		lus_gdp=log(gdp[param$fips=="US"])  #re-initialize US gdp
		lgdp=log(gdp)
		poor=(gdp<3170)  #redefine who's poor
		}
	mgdp=mean(gdp[inipoor])
	if (per[tm]=="mid") {
		impact[b,tm]=(mgdp-mgdp_nocc_mid)/mgdp_nocc_mid*100
		} else {
		impact[b,tm]=(mgdp-mgdp_nocc_end)/mgdp_nocc_end*100
		}
}
print(b)
}

dimnames(impact)[[2]]=names(temp)[3:104]
save(impact,file="DJO_projections.Rdata")
